#setwd('/afs/inf.ed.ac.uk/user/s17/s1725186/Documents/PhD-Models/FirstPUModel/RMarkdowns')

library(tidyverse) ; library(reshape2) ; library(glue) ; library(plotly) ; library(plotlyutils)
library(RColorBrewer) ; library(viridis) ; require(gridExtra) ; library(GGally)
library(Rtsne)
library(ClusterR)
library(DESeq2)

Load preprocessed dataset (preprocessing code in 19_10_14_data_preprocessing.Rmd)

# Gandal dataset
load('./../Data/Gandal/preprocessed_data.RData')
datExpr = datExpr %>% data.frame
DE_info = DE_info %>% data.frame

# GO Neuronal annotations
GO_annotations = read.csv('./../../../PhD-InitialExperiments/FirstYearReview/Data/GO_annotations/genes_GO_annotations.csv')
GO_neuronal = GO_annotations %>% filter(grepl('neuron', go_term)) %>% 
              mutate('ID'=as.character(ensembl_gene_id)) %>% 
              dplyr::select(-ensembl_gene_id) %>% distinct(ID) %>%
              mutate('Neuronal'=1)

# SFARI Genes
SFARI_genes = read_csv('./../../../PhD-InitialExperiments/FirstYearReview/Data/SFARI/SFARI_genes_08-29-2019_with_ensembl_IDs.csv')

# Update DE_info with SFARI and Neuronal information
DE_info = DE_info %>% mutate('ID'=rownames(.)) %>% left_join(SFARI_genes, by='ID') %>% 
  mutate(`gene-score`=ifelse(is.na(`gene-score`), 'None', `gene-score`)) %>%
  distinct(ID, .keep_all = TRUE) %>% left_join(GO_neuronal, by='ID') %>%
  mutate(Neuronal=ifelse(is.na(Neuronal), 0, Neuronal)) %>%
  mutate(gene.score=ifelse(`gene-score`=='None' & Neuronal==1, 'Neuronal', `gene-score`), significant=padj<0.05 & !is.na(padj))


All samples together

plot_data = data.frame('ID'=rownames(datExpr), 'Mean'=rowMeans(datExpr))
p1 = ggplotly(plot_data %>% ggplot(aes(Mean)) + geom_density(color='#0099cc', fill='#0099cc', alpha=0.3) +
              scale_x_log10() + theme_minimal())

plot_data = data.frame('ID'=colnames(datExpr), 'Mean'=colMeans(datExpr))
p2 = ggplotly(plot_data %>% ggplot(aes(Mean)) + geom_density(color='#0099cc', fill='#0099cc', alpha=0.3) +
              theme_minimal() + ggtitle('Mean expression density by gene (left) and by Sample (right)'))

subplot(p1, p2, nrows=1)
rm(p1, p2, plot_data)
plot_data = data.frame('ID'=rownames(datExpr), 'Mean'=rowMeans(datExpr)) %>% 
            left_join(GO_neuronal, by='ID') %>% mutate('Neuronal'=ifelse(is.na(Neuronal),F,T))
p1 = plot_data %>% ggplot(aes(Mean, color=Neuronal, fill=Neuronal)) + geom_density(alpha=0.3) +
                   scale_x_log10() + theme_minimal() + theme(legend.position='bottom') + 
                   ggtitle('Mean expression density by gene')

plot_data = data.frame('ID'=colnames(datExpr), 'Mean'=colMeans(datExpr)) %>% 
            mutate('ID'=substring(ID,2)) %>% left_join(datMeta, by=c('ID'='Dissected_Sample_ID'))
p2 = plot_data %>% ggplot(aes(Mean, color=Diagnosis_, fill=Diagnosis_)) + geom_density(alpha=0.3) +
                   theme_minimal() + theme(legend.position='bottom') + 
                   ggtitle('Mean expression density by Sample')

grid.arrange(p1, p2, nrow=1)

rm(GO_annotations, plot_data, p1, p2)


ASD vs CTL

In general there doesn’t seem to be a lot of variance in mean expression between autism and control samples by gene.

plot_data = data.frame('ID'=rownames(datExpr),
                       'ASD'=rowMeans(datExpr[,datMeta$Diagnosis_=='ASD']),
                       'CTL'=rowMeans(datExpr[,datMeta$Diagnosis_!='ASD']))

plot_data %>% ggplot(aes(ASD,CTL)) + geom_point(alpha=0.1, color='#0099cc') + 
         geom_abline(color='gray') + ggtitle('Mean expression ASD vs CTL') + theme_minimal()

plot_data = rbind(data.frame('Mean'=rowMeans(datExpr[,datMeta$Diagnosis_=='ASD']), 'Diagnosis'='ASD'),
                  data.frame('Mean'=rowMeans(datExpr[,datMeta$Diagnosis_!='ASD']), 'Diagnosis'='CTL')) %>%
            mutate('Diagnosis'=factor(Diagnosis, levels=c('CTL','ASD')))
p1 = ggplotly(plot_data %>% ggplot(aes(Mean, color=Diagnosis, fill=Diagnosis)) + 
              geom_density(alpha=0.3) + scale_x_log10() + theme_minimal())

plot_data = rbind(data.frame('Mean'=colMeans(datExpr[,datMeta$Diagnosis_=='ASD']), 'Diagnosis'='ASD'),
                  data.frame('Mean'=colMeans(datExpr[,datMeta$Diagnosis_!='ASD']), 'Diagnosis'='CTL')) %>%
            mutate('Diagnosis'=factor(Diagnosis, levels=c('CTL','ASD')))
p2 = ggplotly(plot_data %>% ggplot(aes(Mean, color=Diagnosis, fill=Diagnosis)) + 
              geom_density(alpha=0.3) + theme_minimal() +
              ggtitle('Mean expression by Gene (left) and by Sample (right) grouped by Diagnosis'))

subplot(p1, p2, nrows=1)
rm(p1, p2, plot_data)




Visualisations


Samples

PCA

The first principal component seems to separate almost perfectly the two diagnosis

pca = datExpr %>% t %>% prcomp

plot_data = data.frame('ID'=colnames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>% 
            mutate('ID'=substring(ID,2)) %>% left_join(datMeta, by=c('ID'='Dissected_Sample_ID')) %>% 
            dplyr::select('ID','PC1','PC2','Diagnosis_') %>% 
            mutate('Diagnosis'=factor(Diagnosis_, levels=c('CTL','ASD')))

plot_data %>% ggplot(aes(PC1, PC2, color=Diagnosis)) + geom_point() + theme_minimal() + ggtitle('PCA') +
              xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1],1),'%)')) +
              ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2],1),'%)'))

rm(pca, plot_data)


MDS

Looks exactly the same as the PCA visualisation, just inverting both axis

d = datExpr %>% t %>% dist
fit = cmdscale(d, k=2)

plot_data = data.frame('ID'=colnames(datExpr), 'C1'=fit[,1], 'C2'=fit[,2]) %>%
            mutate('ID'=substring(ID,2)) %>% left_join(datMeta, by=c('ID'='Dissected_Sample_ID')) %>% 
            dplyr::select('C1','C2','Diagnosis_') %>%
            mutate('Diagnosis'=factor(Diagnosis_, levels=c('CTL','ASD')))

plot_data %>% ggplot(aes(C1, C2, color=Diagnosis)) + geom_point() + theme_minimal() + ggtitle('MDS')

rm(d, fit, plot_data)


t-SNE

Higher perplexities seem to capture the difference between diagnosis better, although at the end they seem to be capturing another pattern as well, since the samples seem to be grouped in pairs

perplexities = c(2,5,10,25)
ps = list()

for(i in 1:length(perplexities)){
  set.seed(123)
  tsne = datExpr %>% t %>% Rtsne(perplexity=perplexities[i])
  plot_data = data.frame('ID'=colnames(datExpr), 'C1'=tsne$Y[,1], 'C2'=tsne$Y[,2]) %>%
              mutate('ID'=substring(ID,2)) %>% left_join(datMeta, by=c('ID'='Dissected_Sample_ID')) %>%
              dplyr::select('C1','C2','Diagnosis_','Subject_ID') %>%
              mutate('Diagnosis'=factor(Diagnosis_, levels=c('CTL','ASD')))
  ps[[i]] = plot_data %>% ggplot(aes(C1, C2, color=Diagnosis)) + geom_point() + theme_minimal() +
            ggtitle(paste0('Perplexity=',perplexities[i])) + theme(legend.position='none')
}

grid.arrange(grobs=ps, nrow=2)

rm(ps, perplexities, tsne, i)


The second pattern t-SNE seems to be capturing is the subject the samples belong to!

ggplotly(plot_data %>% ggplot(aes(C1, C2, color=Subject_ID)) + geom_point(aes(id=Subject_ID)) + theme_minimal() + 
         theme(legend.position='none') + ggtitle('t-SNE Perplexity=20 coloured by Subject ID'))
rm(plot_data)

Genes

PCA

  • First Principal Component explains over 99% of the total variance

  • There’s a really strong correlation between the mean expression of a gene and the 1st principal component

pca = datExpr %>% prcomp

plot_data = data.frame( 'PC1' = pca$x[,1], 'PC2' = pca$x[,2], 'MeanExpr'=rowMeans(datExpr))

plot_data %>% ggplot(aes(PC1, PC2, color=MeanExpr)) + geom_point(alpha=0.3) + theme_minimal() + 
              scale_color_viridis() + ggtitle('PCA') +
              xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1],1),'%)')) +
              ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2],1),'%)'))

rm(pca, plot_data)


MDS

Distance matrix is too heavy to calculate and the resulting distance object is to big to even load (3.4GB)

t-SNE

Higher perplexities seem to capture the difference between diagnosis better, although at the end they seem to be capturing another pattern as well, since the samples seem to be grouped in pairs

perplexities = c(1,2,5,10,50,100)
ps = list()

for(i in 1:length(perplexities)){
  tsne = read.csv(paste0('./../Visualisations/tsne_perplexity_',perplexities[i],'.csv'))
  plot_data = data.frame('C1'=tsne[,1], 'C2'=tsne[,2], 'MeanExpr'=rowMeans(datExpr))
  ps[[i]] = plot_data %>% ggplot(aes(C1, C2, color=MeanExpr)) + geom_point(alpha=0.5) + theme_minimal() +
            scale_color_viridis() + ggtitle(paste0('Perplexity=',perplexities[i])) + theme(legend.position='none')
}

grid.arrange(grobs=ps, nrow=2)

rm(perplexities, ps, tsne, i)


Differential Expression Analysis

  • 4309 genes (~26%) are significant using a threshold of 0.05 for the adjusted p-value and a without a log Fold Change threshold (keeping the null hypothesis \(H_0: lfc=0\))

  • All genes don’t have an adjusted p-value (no NAs)

table(DE_info$padj<0.05, useNA='ifany')
## 
## FALSE  TRUE 
## 12291  4309

There are two genes with really high lfc: ENSG00000188257 (PLA2G2A) and ENSG00000134028 (ADAMDEC1)

p = DE_info %>% ggplot(aes(log2FoldChange, padj, color=significant)) + geom_point(alpha=0.2) + 
    scale_y_sqrt() + xlab('log2 Fold Change') + ylab('Adjusted p-value') + theme_minimal()
ggExtra::ggMarginal(p, type = 'density', color='gray', fill='gray', size=10)

rm(p)
  • There is a clear negative relation between lfc and mean expression, for both differentially expressed and not differentially expressed groups of genes
plot_data = data.frame('ID'=rownames(datExpr), 'meanExpr'=rowMeans(datExpr)) %>% left_join(DE_info, by='ID') %>%
            mutate('statisticallySignificant'=padj<0.05 & !is.na(padj))

plot_data %>% ggplot(aes(meanExpr, abs(log2FoldChange), color=statisticallySignificant)) + geom_point(alpha=0.1) +
              geom_smooth(method='lm') + theme_minimal() + scale_y_sqrt() + theme(legend.position = 'bottom') +
              xlab('Mean Expression') + ylab('abs(lfc)') + ggtitle('Log fold change by level of expression')

  • When filtering for differential expression, the points seem to separate ino two clouds depending on whether they are over or underexpressed

  • The top cloud corresponds to the over expressed genes and the bottom to the under expressed ones

datExpr_DE = datExpr[DE_info$significant,]

pca = datExpr_DE %>% prcomp

plot_data = cbind(data.frame('PC1'=pca$x[,1], 'PC2'=pca$x[,2]), DE_info[DE_info$significant==TRUE,])

pos_zero = -min(plot_data$log2FoldChange)/(max(plot_data$log2FoldChange)-min(plot_data$log2FoldChange))
p = plot_data %>% ggplot(aes(PC1, PC2, color=log2FoldChange)) + geom_point(alpha=0.5) +
                  scale_color_gradientn(colours=c('#F8766D','#faa49e','white','#00BFC4','#009499'), 
                                        values=c(0, pos_zero-0.01, pos_zero, pos_zero+0.01, 1)) +
                  theme_minimal() + ggtitle('
PCA of differentially expressed genes') + # This is on purpose, PDF doesn't save well without this white space (?)
                  xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1],1),'%)')) +
                  ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2],1),'%)')) + theme(legend.position = 'bottom')
ggExtra::ggMarginal(p, type='density', color='gray', fill='gray', size=10)

rm(pos_zero, p)

Separating the genes into these two groups: Salmon: under-expressed, aqua: over-expressed

  • I separate them with a line instead of by group because in previous preprocessings the relation between the clouds and the over/underexpression of the genes wasn’t as clear as here
intercept=-0.1
slope=0.025

plot_data = plot_data %>% mutate('Group'=ifelse(PC2>slope*PC1+intercept,'overexpressed','underexpressed')) %>%
            mutate('Group' = factor(Group, levels=c('underexpressed','overexpressed')))
plot_data %>% ggplot(aes(PC1, PC2, color=Group)) + geom_point(alpha=0.4) + 
              geom_abline(aes(intercept=intercept, slope=slope), color='gray') +
              xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1],1),'%)')) +
              ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2],1),'%)')) +
              theme_minimal() + ggtitle('PCA')

rm(pca)

Plotting the mean expression by group they seem to have two and three underlying distributions, respectively, so a Gaussian Mixture Model was fitted to each one to separate them into two/three Gaussians and then the points corresponding to each one were plotted in the original PCA plot.

  • Under-expressed ASD genes tend to have a higher mean expression than over-expressed ones (I would have thought this would be the opposite way)
gg_colour_hue = function(n) {
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
}

tot_n_clusters = 5

plot_data = plot_data %>% mutate('MeanExpr'=rowMeans(datExpr_DE), 'SDExpr'=apply(datExpr_DE,1,sd))

GMM_G1 = plot_data %>% filter(Group=='overexpressed') %>% dplyr::select(MeanExpr) %>% GMM(3)
GMM_G2 = plot_data %>% filter(Group=='underexpressed') %>% dplyr::select(MeanExpr) %>% GMM(2)

memberships_G1 = data.frame('ID'=plot_data$ID[plot_data$Group=='overexpressed'],
                            'Membership'=GMM_G1$Log_likelihood %>% apply(1, function(x) glue('over_', which.max(x))))
memberships_G2 = data.frame('ID'=plot_data$ID[plot_data$Group=='underexpressed'],
                            'Membership'=GMM_G2$Log_likelihood %>% apply(1, function(x) glue('under_', which.max(x))))

plot_data = rbind(memberships_G1, memberships_G2) %>% left_join(plot_data, by='ID')
## Warning: Column `ID` joining factor and character vector, coercing into
## character vector
p1 = plot_data %>% ggplot(aes(x=MeanExpr, color=Group, fill=Group)) + geom_density(alpha=0.4) + 
     theme_minimal() + theme(legend.position='bottom')

p2 = plot_data %>% ggplot(aes(x=MeanExpr)) +
     stat_function(fun=dnorm, n=100, colour=gg_colour_hue(tot_n_clusters)[1],
                   args=list(mean=GMM_G1$centroids[1], sd=GMM_G1$covariance_matrices[1])) +
     stat_function(fun=dnorm, n=100, colour=gg_colour_hue(tot_n_clusters)[2],
                   args=list(mean=GMM_G1$centroids[2], sd=GMM_G1$covariance_matrices[2])) +
     stat_function(fun=dnorm, n=100, colour=gg_colour_hue(tot_n_clusters)[3],
                   args=list(mean=GMM_G1$centroids[3], sd=GMM_G1$covariance_matrices[3])) +
     stat_function(fun=dnorm, n=100, colour=gg_colour_hue(tot_n_clusters)[4],
                   args=list(mean=GMM_G2$centroids[1], sd=GMM_G2$covariance_matrices[1])) +
     stat_function(fun=dnorm, n=100, colour=gg_colour_hue(tot_n_clusters)[5],
                  args=list(mean=GMM_G2$centroids[2], sd=GMM_G2$covariance_matrices[2])) +
     theme_minimal()

p3 = plot_data %>% ggplot(aes(PC1, PC2, color=Membership)) + geom_point(alpha=0.4) + theme_minimal() + 
     theme(legend.position='bottom')

grid.arrange(p1, p2, p3, nrow=1)

rm(gg_color_hue, n_clusters, GMM_G1, GMM_G2, memberships_G1, memberships_G2, p1, p2, p3, tot_n_clusters)
## Warning in rm(gg_color_hue, n_clusters, GMM_G1, GMM_G2, memberships_G1, :
## object 'gg_color_hue' not found
## Warning in rm(gg_color_hue, n_clusters, GMM_G1, GMM_G2, memberships_G1, :
## object 'n_clusters' not found

For previous preprocessing pipelines, the pattern found above was also present in the standard deviation, but there doesn’t seem to be any distinguishable patterns now. This could be because the variance was almost homogenised with the vst normalisation algorithm.

plot_data %>% ggplot(aes(x=SDExpr, color=Group, fill=Group)) + geom_density(alpha=0.4) + theme_minimal()

rm(plot_data)



Effects of modifying the log fold change threshold

fc_list = seq(1, 1.2, 0.01)

n_genes = nrow(datExpr)

# Calculate PCAs
datExpr_pca_samps = datExpr %>% data.frame %>% t %>% prcomp(scale.=TRUE)
datExpr_pca_genes = datExpr %>% data.frame %>% prcomp(scale.=TRUE)

# Initialice DF to save PCA outputs
pcas_samps = datExpr_pca_samps$x %>% data.frame %>% dplyr::select(PC1:PC2) %>% 
             mutate('ID'=colnames(datExpr), 'fc'=0, PC1=scale(PC1), PC2=scale(PC2))
pcas_genes = datExpr_pca_genes$x %>% data.frame %>% dplyr::select(PC1:PC2) %>% 
             mutate('ID'=rownames(datExpr), 'fc'=0, PC1=scale(PC1), PC2=scale(PC2))

pca_samps_old = pcas_samps
pca_genes_old = pcas_genes

for(fc in fc_list){
  
  # Recalculate DE_info with the new threshold (p-values change) an filter DE genes
  DE_genes = results(dds, lfcThreshold=log2(fc), altHypothesis='greaterAbs') %>% data.frame %>%
             mutate('ID'=rownames(.)) %>% filter(padj<0.05)
  
  datExpr_DE = datExpr %>% data.frame %>% filter(rownames(.) %in% DE_genes$ID)
  n_genes = c(n_genes, nrow(DE_genes))
  
  # Calculate PCAs
  datExpr_pca_samps = datExpr_DE %>% t %>% prcomp(scale.=TRUE)
  datExpr_pca_genes = datExpr_DE %>% prcomp(scale.=TRUE)

  # Create new DF entries
  pca_samps_new = datExpr_pca_samps$x %>% data.frame %>% dplyr::select(PC1:PC2) %>% 
                  mutate('ID'=colnames(datExpr), 'fc'=fc, PC1=scale(PC1), PC2=scale(PC2))
  pca_genes_new = datExpr_pca_genes$x %>% data.frame %>% dplyr::select(PC1:PC2) %>% 
                  mutate('ID'=DE_genes$ID, 'fc'=fc, PC1=scale(PC1), PC2=scale(PC2))  
  
  # Change PC sign if necessary
  if(cor(pca_samps_new$PC1, pca_samps_old$PC1)<0) pca_samps_new$PC1 = -pca_samps_new$PC1
  if(cor(pca_samps_new$PC2, pca_samps_old$PC2)<0) pca_samps_new$PC2 = -pca_samps_new$PC2
  if(cor(pca_genes_new[pca_genes_new$ID %in% pca_genes_old$ID,]$PC1,
         pca_genes_old[pca_genes_old$ID %in% pca_genes_new$ID,]$PC1)<0){
    pca_genes_new$PC1 = -pca_genes_new$PC1
  }
  if(cor(pca_genes_new[pca_genes_new$ID %in% pca_genes_old$ID,]$PC2, 
         pca_genes_old[pca_genes_old$ID %in% pca_genes_new$ID,]$PC2 )<0){
    pca_genes_new$PC2 = -pca_genes_new$PC2
  }
  
  pca_samps_old = pca_samps_new
  pca_genes_old = pca_genes_new
  
  # Update DFs
  pcas_samps = rbind(pcas_samps, pca_samps_new)
  pcas_genes = rbind(pcas_genes, pca_genes_new)
  
}

# Add Diagnosis/SFARI score information
pcas_samps = pcas_samps %>% mutate('ID'=substring(ID,2)) %>% 
             left_join(datMeta, by=c('ID'='Dissected_Sample_ID')) %>%
             dplyr::select(ID, PC1, PC2, fc, Diagnosis_, Brain_lobe)
# pcas_genes = pcas_genes %>% left_join(SFARI_genes, by='ID') %>% 
#              mutate('score'=as.factor(`gene-score`)) %>%
#              dplyr::select(ID, PC1, PC2, lfc, score)

# Plot change of number of genes
ggplotly(data.frame('fc'=fc_list, 'n_genes'=n_genes[-1]) %>% ggplot(aes(x=fc, y=n_genes)) + 
         geom_point() + geom_line() + theme_minimal() + xlab('Fold Change') +
         ggtitle('Number of remaining genes when modifying filtering threshold'))
rm(fc_list, n_genes, fc, pca_samps_new, pca_genes_new, pca_samps_old, pca_genes_old, 
   datExpr_pca_samps, datExpr_pca_genes)


Samples

Note: PC values get smaller as Log2 fold change increases, so on each iteration the values were scaled so it would be easier to compare between frames

Coloured by Diagnosis:

Fold changes between 1 and 1.1 seem to separate the two diagnosis groups best

ggplotly(pcas_samps %>% ggplot(aes(PC1, PC2, color=Diagnosis_)) + geom_point(aes(frame=fc, ids=ID)) + 
         theme_minimal() + ggtitle('Samples PCA plot modifying filtering threshold'))


Coloured by brain region:

There doesn’t seem to be any recognisable pattern

ggplotly(pcas_samps %>% ggplot(aes(PC1, PC2, color=Brain_lobe)) + geom_point(aes(frame=fc, ids=ID)) + 
         theme_minimal() + ggtitle('Samples PCA plot modifying filtering threshold'))


Genes

ggplotly(pcas_genes %>% ggplot(aes(PC1, PC2)) + geom_point(aes(frame=fc, ids=ID, alpha=0.3)) + 
         theme_minimal() + ggtitle('Genes PCA plot modifying filtering threshold'))




Session info

sessionInfo()
## R version 3.5.2 (2018-12-20)
## Platform: x86_64-redhat-linux-gnu (64-bit)
## Running under: Scientific Linux 7.6 (Nitrogen)
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
## 
## locale:
##  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] DESeq2_1.22.2               SummarizedExperiment_1.12.0
##  [3] DelayedArray_0.8.0          BiocParallel_1.16.6        
##  [5] matrixStats_0.54.0          Biobase_2.42.0             
##  [7] GenomicRanges_1.34.0        GenomeInfoDb_1.18.2        
##  [9] IRanges_2.16.0              S4Vectors_0.20.1           
## [11] BiocGenerics_0.28.0         ClusterR_1.1.8             
## [13] gtools_3.5.0                Rtsne_0.15                 
## [15] GGally_1.4.0                gridExtra_2.3              
## [17] viridis_0.5.1               viridisLite_0.3.0          
## [19] RColorBrewer_1.1-2          plotlyutils_0.0.0.9000     
## [21] plotly_4.8.0                glue_1.3.1                 
## [23] reshape2_1.4.3              forcats_0.3.0              
## [25] stringr_1.4.0               dplyr_0.8.0.1              
## [27] purrr_0.3.1                 readr_1.3.1                
## [29] tidyr_0.8.3                 tibble_2.1.1               
## [31] ggplot2_3.1.0               tidyverse_1.2.1            
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.1.0           backports_1.1.2        Hmisc_4.1-1           
##   [4] plyr_1.8.4             lazyeval_0.2.2         splines_3.5.2         
##   [7] gmp_0.5-13.5           crosstalk_1.0.0        digest_0.6.18         
##  [10] htmltools_0.3.6        magrittr_1.5           checkmate_1.8.5       
##  [13] memoise_1.1.0          cluster_2.0.7-1        annotate_1.60.1       
##  [16] modelr_0.1.4           colorspace_1.4-1       blob_1.1.1            
##  [19] rvest_0.3.2            haven_1.1.1            xfun_0.5              
##  [22] crayon_1.3.4           RCurl_1.95-4.10        jsonlite_1.6          
##  [25] genefilter_1.64.0      survival_2.43-3        ape_5.3               
##  [28] gtable_0.2.0           zlibbioc_1.28.0        XVector_0.22.0        
##  [31] abind_1.4-5            scales_1.0.0           DBI_1.0.0             
##  [34] miniUI_0.1.1.1         Rcpp_1.0.1             xtable_1.8-3          
##  [37] htmlTable_1.11.2       magic_1.5-9            foreign_0.8-71        
##  [40] bit_1.1-14             Formula_1.2-3          htmlwidgets_1.3       
##  [43] httr_1.4.0             acepack_1.4.1          pkgconfig_2.0.2       
##  [46] reshape_0.8.7          XML_3.98-1.11          nnet_7.3-12           
##  [49] locfit_1.5-9.1         tidyselect_0.2.5       labeling_0.3          
##  [52] rlang_0.3.2            later_0.8.0            AnnotationDbi_1.44.0  
##  [55] munsell_0.5.0          cellranger_1.1.0       tools_3.5.2           
##  [58] cli_1.1.0              generics_0.0.2         RSQLite_2.1.1         
##  [61] ade4_1.7-11            broom_0.5.1            evaluate_0.13         
##  [64] geometry_0.4.0         FD_1.0-12              yaml_2.2.0            
##  [67] knitr_1.22             bit64_0.9-7            nlme_3.1-137          
##  [70] mime_0.6               ggExtra_0.8            xml2_1.2.0            
##  [73] compiler_3.5.2         rstudioapi_0.7         geneplotter_1.60.0    
##  [76] stringi_1.4.3          lattice_0.20-38        Matrix_1.2-15         
##  [79] vegan_2.5-2            permute_0.9-4          pillar_1.3.1          
##  [82] data.table_1.12.0      bitops_1.0-6           httpuv_1.5.0          
##  [85] R6_2.4.0               latticeExtra_0.6-28    promises_1.0.1        
##  [88] MASS_7.3-51.1          assertthat_0.2.1       withr_2.1.2           
##  [91] GenomeInfoDbData_1.2.0 mgcv_1.8-26            hms_0.4.2             
##  [94] grid_3.5.2             rpart_4.1-13           rmarkdown_1.12        
##  [97] Cairo_1.5-9            shiny_1.2.0            lubridate_1.7.4       
## [100] base64enc_0.1-3